***************************************************************************************************************************
/*	This do-file
			a. The distribution of impact standard deviations for testing the null og zero impact standard deviation as in Heckman & Smith & Clements, 1999. 
				1. Sample (ssize) from control distribution reps times. 
				2. Sample (ssize) from control distribution reps times and add mtimpact to create full treatment and cctimpact to create reduced treatment. 
				3. Collapses all the distributions into percentiles. 
				4. Compute differences between treatments and controls for each percentile. 
				5. Compute the standard deviation of the percentile differences.    
				6. The command fhbounds is define with the syntax : syntax varlist(max=1) [if], Reps(# of repetitions) MTimpact(assumed value of mt impact) CCTimpact(assumed value of mt impact) ssize(size of the samples)
				
				
	REMEMBER TO DETAIL AND EXPLAIN WHAT THE DO-FILE DOES!!!			
*/
***************************************************************************************************************************
cap program drop isd_pval
program define isd_pval
syntax varlist(max=1) [if], Reps(integer) MTimpact(real) CCTimpact(real) ssize_c(integer) ssize_t(integer)  [se BSReps(integer 50)]
cap mat drop store 
scalar drop _all
	// Create Storing Matrices
	if "`se'"!="" {
		display "Number of boostrap samples assumed: `bsreps'"
		local se_size = 8
	}
	else {
		local se_size = 2
	}
	mat store = J(1,`se_size',.)
	// Temporary save Data
	tempfile data
	qui save `data'
	display "MT assumed impact = `mtimpact'"
	display "CCT assumed impact = `cctimpact'"
	// Start Monte Carlo Simulations
	local counter = 0
	set seed 93478	
	while `counter' < `reps' { 
	use `data', clear
	// Keep controls
		qui keep if Study_Arm == 0
		if "`if'" != "" {
			keep `if'
		}
		keep `varlist'
		// Temporary save Control
		tempfile controls
		qui save `controls'
	// Sample from control Group
		use `controls', clear
		bsample `ssize_c'
		// Collapse into percentiles
		pctile p_control  = `varlist', nq(100)
		keep p_control
		qui keep if p_control!=.
		sort p_control
		gen p = _n
		tempfile control_sample
		qui save `control_sample'
	// Sample MT and CCT separately: form treated group by adding X and sampling from control group
		use `controls', clear
		foreach t in mt cct {
			use `controls', clear
			qui sum `varlist'
			local mean1 = r(mean)
			expand 10
			qui sum `varlist'
			local mean2 = r(mean)
			// Check that mean is not modified
			assert `mean1'==`mean2'
			// Sample with replacement
			bsample `ssize_t'
			gen `t'_impact  = `varlist' + ``t'impact'
			// Collapse into percentiles
			pctile p_`t'  = `t'_impact, nq(100)
			keep p_`t'
			qui keep if p_`t'!=. 
			sort p_`t'
			gen p = _n
			// Save
			tempfile `t'_sample
			qui save ``t'_sample'
		}	
	// Join controls with treatments
		use `control_sample', clear
		qui merge 1:1 p using `mt_sample'
		drop _merge
		qui merge 1:1 p using `cct_sample'
		drop _merge
	// Take percentile differences:  "The standard deviation of the differences is the estimated impact standard deviation for a given sample"
		gen p_diff_mt  = p_mt - p_control
		gen p_diff_cct = p_cct - p_control
	// SD of percentile differences
		qui sum p_diff_mt 
		scalar ISD_MT = r(sd) 
		qui sum p_diff_cct 
		scalar ISD_CCT = r(sd)
	// Boostrap Standard error: "The standard deviation of the estimated impact standard deviations in the bootstrap samples is the bootstrap standard error for the corresponding sample"
		if "`se'"!="" {
			qui bootstrap r(sd),  reps(`bsreps'): sum p_diff_mt 
			mat m = e(se)
			scalar ISD_MT_se = m[1,1] 
			qui bootstrap r(sd),  reps(`bsreps'): sum p_diff_cct
			mat c = e(se)
			scalar ISD_CCT_se = c[1,1] 
			// Bootstrap Confidence interval:  Centred on the estimated impact standard deviation for each data sample and extending 1.96 times the bootstraph s.e. on each side
			scalar ubs_ci_MT  = ISD_MT + (1.96*ISD_MT_se)
			scalar lbs_ci_MT  = ISD_MT - (1.96*ISD_MT_se)
			scalar ubs_ci_CCT = ISD_CCT + (1.96*ISD_CCT_se)
			scalar lbs_ci_CCT = ISD_CCT - (1.96*ISD_CCT_se)
		}
		// Display replication
		local j = `counter'+1
		display "ISD in replication `j':"
		scalar list ISD_MT ISD_CCT
		local counter = `counter'+1
		local j = `j'+1
	// Store values
		if "`se'"!="" {
			mat store = store\(ISD_MT,ISD_CCT,ISD_MT_se,ISD_CCT_se,ubs_ci_MT,lbs_ci_MT,ubs_ci_CCT,lbs_ci_CCT)
		}
		else {
			mat store = store\(ISD_MT,ISD_CCT)
		}
	}	
	mat store = store[2..`j',1..`se_size']
	mat list store
	// Compute the Distribution 
	clear
	svmat store
	ren store1 isd_MT
	ren store2 isd_CCT
	cap ren store3 isd_MT_se
	cap ren store4 isd_CCT_se
	cap ren store5 ubs_ci_MT 
	cap ren store6 lbs_ci_MT 
	cap ren store7 ubs_ci_CCT 
	cap ren store8 lbs_ci_CCT
	// Collapse into new distributions by percentiles
end